LTER Pulses WaveletComp R Package Tutorial
KM Hall
2022-09-16
The following is a tutorial for the 2022 LTER All Scientists’ Meeting pulse dynamics workshop about the basic use of the WaveletComp R package for wavelet and cross-wavelet analysis of time series data of one and two variables, respectively.
For a thorough discussion of WaveletComp, refer to the excellent guided tour document that accompanies the package. The guided tour goes into much greater detail than this tutorial, and much of this tutorial’s materials are based on the guided tour.
A simple, template R script is also available along with this tutorial - see LTER_Pulses_WaveletComp_template.R.
Basic background on wavelets
Wavelets provide a method for studying periodic time series phenomena that may be either stationary or non-stationary over time. This allows for the examination of time series data in both the time and frequency domains.
A wavelet function, the “mother” wavelet, is scaled (stretched and compressed) into a series of wavelet children. The children are shifted/translated across the entire time series. This is called a convolution of the time series with the wavelet children and is the wavelet transform.
There are many types of wavelets and WaveletComp uses the Morlet wavelet:
The wavelet analysis provides the wavelet power spectrum that can be interpreted as time-frequency wavelet energy density, and is related to the amplitude of a periodic signal at a given point in time.
The relationship between period and frequency is:
Period = 1 / Frequency
There is a trade off between time and frequency resolution. Long periods are well detected in the frequency domain, but poorly localized in the time domain. Short periods are well detected in the time domain, but poorly localized in the frequency domain.
Through the wavelet transform, it is possible to identify what periodicities are most strongly present in a univariate time series.
By utilizing a series of simulations, WaveletComp tests the null hypothesis of no periodicity using a white noise process or several other test options. Regions of the time series that are significantly different from the null may be identified in the model output and visualizations.
Cross-wavelet analysis allows for comparison of the frequency contents of a bivariate time series, thus allowing for an examination of synchronicity or coherency between the two time series. Coherency is analogous to classical correlation. The cross-wavelet power shows the similarity of the bivariate time series’ wavelet power in the time-frequency domain. It also allows for the assessment of how closely in/out of phase the two time series are.
Some examples will help to provide a better understanding of the wavelet transform of time series data.
The Data
One synthetic data example is provided that is directly from the guided tour.
The other data used in the examples below are daily AmeriFlux data from the US-Seg: Sevilleta grassland flux tower located on the Sevilleta National Wildlife Refuge in New Mexico from 2007-2020. Latitude/Longitude: 34.3623, -106.7020. Elevation: 1596 m. See the AmeriFlux site info page for more information. Variables from the data that are used in the tutorial include daily mean air temperature, precipitation, and net ecosystem exchange (NEE).
The Context
Utilize wavelets to conduct cross-site comparisons of pulse dynamics within (and beyond) the LTER Network.
Motivating questions include:
- At what temporal scales are the major sources of variance
(pulsed-ness) in abiotic resources and conditions?
- How have the pulse regimes of abiotic resources and conditions
changed over time?
- How strongly have ecological responses tracked abiotic pulse regimes, and how large are the lags in response?
The Setup
Several R packages are loaded along with the daily US-Seg data.
library(WaveletComp)
library(dplyr)
library(matrixStats)
library(tidyr)
library(ggplot2)
library(lubridate)
library(readr)seg <- read_csv("~/Documents/LTER_Pulse_Time_Series_WG/SEV_daily_flux/US-Seg_daily_aflx.csv") %>%
mutate(date = mdy(Date)) %>% # WaveletComp expects date var named "date"
select(-Date) %>%
select(date, everything()) %>%
arrange(date) # make sure data are properly sorted in chronological order## Rows: 5114 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (22): GPP_int, GPP_dayint, RECO_int, RECO_dayint, RECO_nightint, NEE_int...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Let’s take a quick look at the data we are working with.
Data structure:
str(seg)## tibble [5,114 × 23] (S3: tbl_df/tbl/data.frame)
## $ date : Date[1:5114], format: "2007-01-01" "2007-01-02" ...
## $ GPP_int : num [1:5114] 0 0.00395 0.0362 0.25915 0.46783 ...
## $ GPP_dayint : num [1:5114] 0 0.00228 0.0362 0.13995 0.24107 ...
## $ RECO_int : num [1:5114] 0.495 0.61 0.7 0.592 0.514 ...
## $ RECO_dayint : num [1:5114] 0.247 0.262 0.234 0.273 0.268 ...
## $ RECO_nightint: num [1:5114] 0.249 0.348 0.465 0.32 0.246 ...
## $ NEE_int : num [1:5114] 0.495 0.606 0.664 0.333 0.046 ...
## $ NEE_dayint : num [1:5114] 0.2467 0.2594 0.1982 0.1328 0.0272 ...
## $ NEE_nightint : num [1:5114] 0.2486 0.3464 0.4654 0.2005 0.0188 ...
## $ ET_int : num [1:5114] 0.368 0.25 0.255 0.458 0.864 ...
## $ ET_dayint : num [1:5114] 0.359 0.262 0.305 0.425 0.794 ...
## $ P_int : num [1:5114] 0.614 1.024 0 0 1.024 ...
## $ TA_avg : num [1:5114] -8.53 -5.98 -7.56 -1.98 3.52 ...
## $ TA_max : num [1:5114] -2.14 -3.26 -2.93 6.94 9.91 ...
## $ RH_avg : num [1:5114] 88.5 86.5 87.9 72.6 63.1 ...
## $ SW_IN_avg : num [1:5114] 125 128 145 150 150 ...
## $ NETRAD_avg : num [1:5114] -34.28 -33.93 -8.28 -7.74 3.35 ...
## $ PPFD_IN_avg : num [1:5114] 243 249 282 291 291 ...
## $ PPFD_IN_sum : num [1:5114] 11407 11948 13517 13951 13985 ...
## $ VPD_avg : num [1:5114] 0.0375 0.0569 0.0477 0.1822 0.3185 ...
## $ VPD_max : num [1:5114] 0.226 0.267 0.174 0.543 0.774 ...
## $ LE_avg : num [1:5114] 13.28 8.77 9.01 15.69 28.8 ...
## $ H_avg : num [1:5114] 15.94 7 10.08 2.15 -6.45 ...Data summary:
summary(seg)## date GPP_int GPP_dayint RECO_int
## Min. :2007-01-01 Min. :0.0000 Min. :0.00000 Min. :0.001662
## 1st Qu.:2010-07-02 1st Qu.:0.1340 1st Qu.:0.09379 1st Qu.:0.195740
## Median :2013-12-31 Median :0.2662 Median :0.22001 Median :0.345905
## Mean :2013-12-31 Mean :0.5864 Mean :0.53488 Mean :0.521193
## 3rd Qu.:2017-07-01 3rd Qu.:0.7426 3rd Qu.:0.68267 3rd Qu.:0.658740
## Max. :2020-12-31 Max. :5.2217 Max. :5.10939 Max. :5.882377
##
## RECO_dayint RECO_nightint NEE_int NEE_dayint
## Min. :0.0000 Min. :0.00000 Min. :-3.008381 Min. :-3.78218
## 1st Qu.:0.1021 1st Qu.:0.07787 1st Qu.:-0.195228 1st Qu.:-0.34236
## Median :0.1836 Median :0.15992 Median :-0.006409 Median :-0.07767
## Mean :0.2997 Mean :0.22151 Mean :-0.065162 Mean :-0.23520
## 3rd Qu.:0.3883 3rd Qu.:0.27416 3rd Qu.: 0.141722 3rd Qu.: 0.04146
## Max. :4.3728 Max. :1.60844 Max. : 5.810892 Max. : 4.31067
##
## NEE_nightint ET_int ET_dayint P_int
## Min. :-0.85730 Min. :0.01599 Min. :0.0000 Min. : 0.0000
## 1st Qu.: 0.03406 1st Qu.:0.19703 1st Qu.:0.1950 1st Qu.: 0.0000
## Median : 0.11403 Median :0.38021 Median :0.3622 Median : 0.0000
## Mean : 0.17004 Mean :0.58996 Mean :0.5605 Mean : 0.6319
## 3rd Qu.: 0.22711 3rd Qu.:0.75909 3rd Qu.:0.7085 3rd Qu.: 0.0000
## Max. : 1.56843 Max. :4.26213 Max. :4.2112 Max. :46.3441
##
## TA_avg TA_max RH_avg SW_IN_avg
## Min. :-18.914 Min. :-16.44 Min. : 6.062 Min. : 9.236
## 1st Qu.: 6.075 1st Qu.: 13.47 1st Qu.:26.113 1st Qu.:163.145
## Median : 14.267 Median : 21.77 Median :38.656 Median :235.172
## Mean : 13.712 Mean : 20.88 Mean :40.155 Mean :231.808
## 3rd Qu.: 21.993 3rd Qu.: 28.99 3rd Qu.:52.572 3rd Qu.:300.307
## Max. : 30.313 Max. : 38.58 Max. :98.833 Max. :385.824
##
## NETRAD_avg PPFD_IN_avg PPFD_IN_sum VPD_avg
## Min. :-86.76 Min. : 0.0535 Min. : 0 Min. :0.006617
## 1st Qu.: 40.06 1st Qu.: 303.3559 1st Qu.:14138 1st Qu.:0.554238
## Median : 77.19 Median : 442.0772 Median :20836 Median :1.106678
## Mean : 77.84 Mean : 438.4796 Mean :20492 Mean :1.272664
## 3rd Qu.:112.63 3rd Qu.: 576.5788 3rd Qu.:27422 3rd Qu.:1.830674
## Max. :654.11 Max. :1977.1261 Max. :34865 Max. :4.143947
## NA's :195 NA's :100
## VPD_max LE_avg H_avg
## Min. :0.01723 Min. : 0.5317 Min. :-30.02
## 1st Qu.:1.18398 1st Qu.: 6.3996 1st Qu.: 30.39
## Median :2.16125 Median : 12.0348 Median : 50.88
## Mean :2.33556 Mean : 18.5688 Mean : 52.94
## 3rd Qu.:3.34840 3rd Qu.: 23.8515 3rd Qu.: 74.22
## Max. :6.43246 Max. :134.6514 Max. :143.23
## Notice that there are no missing data. WaveletComp requires a complete time series without missing data. It may be necessary for you to preprocess your data to gap fill any missing data.
Data head:
head(seg, 10)## # A tibble: 10 × 23
## date GPP_int GPP_dayint RECO_int RECO_dayint RECO_nightint NEE_int
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2007-01-01 0 0 0.495 0.247 0.249 0.495
## 2 2007-01-02 0.00395 0.00228 0.610 0.262 0.348 0.606
## 3 2007-01-03 0.0362 0.0362 0.700 0.234 0.465 0.664
## 4 2007-01-04 0.259 0.140 0.592 0.273 0.320 0.333
## 5 2007-01-05 0.468 0.241 0.514 0.268 0.246 0.0460
## 6 2007-01-06 0.0969 0.0526 0.374 0.184 0.190 0.277
## 7 2007-01-07 0.115 0.113 0.448 0.189 0.259 0.333
## 8 2007-01-08 0.118 0.104 0.373 0.145 0.228 0.255
## 9 2007-01-09 0.129 0.108 0.314 0.117 0.197 0.185
## 10 2007-01-10 0.196 0.123 0.290 0.146 0.144 0.0938
## # … with 16 more variables: NEE_dayint <dbl>, NEE_nightint <dbl>, ET_int <dbl>,
## # ET_dayint <dbl>, P_int <dbl>, TA_avg <dbl>, TA_max <dbl>, RH_avg <dbl>,
## # SW_IN_avg <dbl>, NETRAD_avg <dbl>, PPFD_IN_avg <dbl>, PPFD_IN_sum <dbl>,
## # VPD_avg <dbl>, VPD_max <dbl>, LE_avg <dbl>, H_avg <dbl>Example 1 - US-Seg Average Daily Temperature Wavelet Analysis
This first example examines the underlying periodic structures of US-Seg mean daily temperature in Celsius. The example exhibits simple periodic structure, but it is a useful introduction to wavelets.
The US-Seg Average Daily Temperature time series:
seg %>%
ggplot(aes(x = date, y = TA_avg)) +
geom_line(size = 0.4, color = "tan") +
labs(title = "US-Seg - Average Daily Temperature\n2007-2020",
x = "Date",
y = "Average Temperature (C)")
The analyze.wavelet function performs the wavelet
transformation and provides the power spectrum for mean daily
temperature for the time series:
seg_ta_avg.w <- analyze.wavelet(seg, "TA_avg",
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE, method = "white.noise",
n.sim = 25) # should set n.sim higher - default is 100## Smoothing the time series...
## Starting wavelet transformation...
## ... and simulations...
##
|
| | 0%
|
|=== | 4%
|
|====== | 8%
|
|======== | 12%
|
|=========== | 16%
|
|============== | 20%
|
|================= | 24%
|
|==================== | 28%
|
|====================== | 32%
|
|========================= | 36%
|
|============================ | 40%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================== | 60%
|
|============================================= | 64%
|
|================================================ | 68%
|
|================================================== | 72%
|
|===================================================== | 76%
|
|======================================================== | 80%
|
|=========================================================== | 84%
|
|============================================================== | 88%
|
|================================================================ | 92%
|
|=================================================================== | 96%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave Phase Ampl Power Power.avg Power.pval Power.avg.pval Ridge Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzanalyze.wavelet has numerous arguments that you can set
and tweak from their defaults, so check out the documentation for the
function using ?analyze.wavelet in the R console.
What do you think the most powerful period(s) will be for the average daily temperature time series?
# you can define your continuous color palette
# mypalette <- "terrain.colors(100)"
wt.image(seg_ta_avg.w,
col.contour = "white",
# color.palette = mypalette,
lwd = 1.5,
show.date = TRUE,
main = "US-Seg Flux Tower Daily Mean Air Temperature (C)",
legend.params = list(width = 1.8, label.digits = 4, lab = "Power"))
The most powerful periods in the time series are around the annual temperature cycle of ~365 days. However, if you look closely, there are a couple of other periods with high power at certain time ranges in the time series.
The average power of the entire time series can be graphed:
# look at the average power of the series
maximum.level = 1.001*max(seg_ta_avg.w$Power.avg) # See WaveletComp Guided Tour for explanation
wt.avg(seg_ta_avg.w, maximum.level = maximum.level,
legend.coords = "bottomright",
main = "Average Power for US-Seg Mean Daily Air Temperature")
Because the wavelet transform preserves information, it is possible to reconstruct the original time series for the time periods whose average power were found to be significantly different from a white noise structure:
reconstruct(seg_ta_avg.w,
show.date = TRUE,
legend.coords = "bottomright",
ylim = c(-40, 20),
main = "Reconstruction of US-Seg Flux Tower \nDaily Mean Airt Temperature",
col = c('tan', 'black'))## Your input object class is 'analyze.wavelet'...
## Your time series 'TA_avg' will be reconstructed...
## Starting the reconstruction process...
## Original (detrended) and reconstructed series are being plotted...
## Class attributes are accessible through following names:
## series rec.waves loess.span lvl only.coi only.sig siglvl only.ridge rnum.used rescale dt dj Period Scale nc nr axis.1 axis.2 date.format date.tzanalyze.wavelet produces several results variables that
can be used for further analysis. After inspecting the wavelet image, it
may be interesting to look more closely at powerful periods in certain
time ranges. The period information output from
analyze.wavelet may be used to pinpoint and examine
periodicities of interest:
seg_ta_avg.w$Period## [1] 2.000000 2.070530 2.143547 2.219139 2.297397 2.378414
## [7] 2.462289 2.549121 2.639016 2.732081 2.828427 2.928171
## [13] 3.031433 3.138336 3.249010 3.363586 3.482202 3.605002
## [19] 3.732132 3.863745 4.000000 4.141060 4.287094 4.438278
## [25] 4.594793 4.756828 4.924578 5.098243 5.278032 5.464161
## [31] 5.656854 5.856343 6.062866 6.276673 6.498019 6.727171
## [37] 6.964405 7.210004 7.464264 7.727491 8.000000 8.282119
## [43] 8.574188 8.876556 9.189587 9.513657 9.849155 10.196485
## [49] 10.556063 10.928322 11.313708 11.712686 12.125733 12.553346
## [55] 12.996038 13.454343 13.928809 14.420007 14.928528 15.454981
## [61] 16.000000 16.564239 17.148375 17.753112 18.379174 19.027314
## [67] 19.698311 20.392970 21.112127 21.856644 22.627417 23.425371
## [73] 24.251465 25.106691 25.992077 26.908685 27.857618 28.840015
## [79] 29.857056 30.909963 32.000000 33.128478 34.296751 35.506223
## [85] 36.758347 38.054628 39.396621 40.785940 42.224253 43.713288
## [91] 45.254834 46.850742 48.502930 50.213382 51.984153 53.817371
## [97] 55.715236 57.680030 59.714111 61.819925 64.000000 66.256955
## [103] 68.593502 71.012446 73.516695 76.109255 78.793242 81.571880
## [109] 84.448506 87.426576 90.509668 93.701485 97.005860 100.426765
## [115] 103.968307 107.634741 111.430472 115.360059 119.428223 123.639850
## [121] 128.000000 132.513910 137.187003 142.024892 147.033389 152.218511
## [127] 157.586485 163.143760 168.897013 174.853153 181.019336 187.402969
## [133] 194.011721 200.853529 207.936613 215.269482 222.860944 230.720118
## [139] 238.856446 247.279700 256.000000 265.027821 274.374006 284.049785
## [145] 294.066779 304.437021 315.172970 326.287521 337.794025 349.706306
## [151] 362.038672 374.805938 388.023441 401.707058 415.873227 430.538965
## [157] 445.721888 461.440237 477.712892 494.559400 512.000000 530.055641
## [163] 548.748013 568.099570 588.133558 608.874043 630.345940 652.575041
## [169] 675.588050 699.412611 724.077344 749.611876 776.046882 803.414116
## [175] 831.746454 861.077929 891.443777 922.880474 955.425783 989.118801
## [181] 1024.000000 1060.111282 1097.496026 1136.199139 1176.267116 1217.748086
## [187] 1260.691879 1305.150082 1351.176101 1398.825223 1448.154688 1499.223753
## [193] 1552.093764 1606.828232 1663.492908Let’s look at the power of the periods around 365. The 365-day period is between records 151 and 152. The following looks at a broad range of periods between 256 and 512 days.
Convert power data to a matrix:
seg_ta_avg.w_power <- as.matrix(seg_ta_avg.w$Power)The dimensions of the Power matrix have the number of periods as columns and the number of records in the time series as rows:
dim(seg_ta_avg.w$Power)## [1] 195 5114From this, it is possible to pull out various results from
analyze.wavelet. For example, here are the results looking
at the periodicities from 256 to 512:
seg_256_to_512_power <- as.matrix(seg_ta_avg.w_power[141:161,])
dimnames(seg_256_to_512_power) <- list(paste0("period_", round(seg_ta_avg.w$Period[141:161], 2)))
seg_256_to_512_power_t <- as_tibble(t(seg_256_to_512_power))
seg_256_to_512_power_df <- tibble(date = seg$date, seg_256_to_512_power_t)
seg_256_to_512_power_df_long <- seg_256_to_512_power_df %>%
pivot_longer(contains("period_"))
seg_256_to_512_power_df_long %>%
ggplot(aes(x = date, y = value, color = name)) +
geom_line() +
facet_wrap(~ name) +
theme(legend.position="none") +
labs(title = "Average daily temperature periods around 365",
x = "Date",
y = "Power")
Looking closely at the wavelet image, there appears to be some high power periods in 2011 between 20 and 50 day periods:
This is the time series for 2011:
seg %>%
filter(year(date) == 2011) %>%
ggplot(aes(x = date, y = TA_avg)) +
geom_line(size = 0.6) +
labs(title = "2011 Daily Average Temperature (C)")
There was a strong freeze in February 2011, and this shows up in the wavelet transform.
To pull out the 2011 data from analyze.wavelet, identify
the 2011 records:
(start_2011 <- ymd("2011-01-01") - ymd(min(seg$date)))## Time difference of 1461 days(end_2011 <- ymd("2011-12-31") - ymd(min(seg$date)))## Time difference of 1825 daysThen look at power bands for periods 20-50 in 2011:
seg_2011_power <- as.matrix(seg_ta_avg.w_power[68:94, 1461:1825])
dimnames(seg_2011_power) <- list(paste0("period_", round(seg_ta_avg.w$Period[68:94], 2)))
seg_2011_power_t <- as_tibble(t(seg_2011_power))
seg_2011_power_df <- tibble(date = seq(from = ymd("2011-01-01"),
to = ymd("2011-12-31"),
by = "1 day"),
seg_2011_power_t)
seg_2011_power_df_long <- seg_2011_power_df %>%
pivot_longer(contains("period_"))
seg_2011_power_df_long %>%
ggplot(aes(x = date, y = value, color = name)) +
geom_line() +
facet_wrap(~ name) +
theme(legend.position="none") +
labs(title = "Average daily temperature in 2011",
x = "Date",
y = "Power") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
This first example was used to show a basic analytic process for analyzing wavelets and how to identify high power regions of interest in a univariate time series.
Example 2 - Non-stationary synthetic data
This synthetic example shows how wavelets can detect changes in periodicity throughout the time series (time-period domain). This is an advantage over other techniques that cannot detect non-stationarity.
This example comes directly from the guided tour.
First, produce synthetic/pretend data:
# made up data for example
x1 <- periodic.series(start.period = 100, length = 400)
x2 <- 1.5*periodic.series(start.period = 50, length = 200)
x <- c(x1, x2, x1) + 0.2*rnorm(1000)
pretend.data <- data.frame(x = x)pretend.data %>%
ggplot(aes(x = seq(from = 1, to = 1000, by=1), y = x)) +
geom_line(size = 0.6, color = "tan") +
labs(title = "Non-Stationary Data Example",
x = "time",
y= "value")
Then conduct the wavelet transform and visualize:
pretend.w <- analyze.wavelet(pretend.data, "x",
method = "white.noise",
loess.span = 0,
dt = 1, dj = 1/250,
lowerPeriod = 32, upperPeriod = 256,
make.pval = TRUE, n.sim = 10)## Starting wavelet transformation...
## ... and simulations...
##
|
| | 0%
|
|======= | 10%
|
|============== | 20%
|
|===================== | 30%
|
|============================ | 40%
|
|=================================== | 50%
|
|========================================== | 60%
|
|================================================= | 70%
|
|======================================================== | 80%
|
|=============================================================== | 90%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave Phase Ampl Power Power.avg Power.pval Power.avg.pval Ridge Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzwt.image(pretend.w, color.key = "interval", n.levels = 250,
legend.params = list(lab = "wavelet power levels"))
Notice how the wavelet image clearly identifies the changes in dominant periods throughout the time series.
Example 3 - Cross-wavelet analysis - US-Seg Precipitation and Net Ecosystem Exchange
Cross-wavelet analysis can be used to examing two time series variables at the same time to identify synchronicities and phase differences between the variables.
The two variables for this analysis are:
- Daily Precipitation (mm) - P_int
- Net Ecosystem Exchange (gC/m^2) - NEE_int
Before conducting the cross-wavelet analysis, let’s take look at the two variables by themselves.
seg %>%
ggplot(aes(x = date, y = P_int)) +
geom_line(size = 0.4, color = "tan") +
labs(title = "US-Seg - Daily Precipitation",
x = "Date",
y = "Precipitation (mm)")
seg %>%
ggplot(aes(x = date, y = NEE_int)) +
geom_line(size = 0.4, color = "tan") +
labs(title = "US-Seg - Net Ecosystem Exchange",
x = "Date",
y = "NEE (gC/m^2)")
And conduct univariate wavelet transforms:
seg_p_int.w <- analyze.wavelet(seg, "P_int",
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE, method = "white.noise",
n.sim = 25)## Smoothing the time series...
## Starting wavelet transformation...
## ... and simulations...
##
|
| | 0%
|
|=== | 4%
|
|====== | 8%
|
|======== | 12%
|
|=========== | 16%
|
|============== | 20%
|
|================= | 24%
|
|==================== | 28%
|
|====================== | 32%
|
|========================= | 36%
|
|============================ | 40%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================== | 60%
|
|============================================= | 64%
|
|================================================ | 68%
|
|================================================== | 72%
|
|===================================================== | 76%
|
|======================================================== | 80%
|
|=========================================================== | 84%
|
|============================================================== | 88%
|
|================================================================ | 92%
|
|=================================================================== | 96%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave Phase Ampl Power Power.avg Power.pval Power.avg.pval Ridge Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzmypalette <- "terrain.colors(100)"
wt.image(seg_p_int.w,
col.contour = "white",
color.palette = mypalette,
lwd = 1,
show.date = TRUE,
main = "US-Seg Flux Tower Daily Precipitation (mm)",
legend.params = list(width = 1.8, label.digits = 4))
seg_nee_int.w <- analyze.wavelet(seg, "NEE_int",
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE, method = "white.noise",
n.sim = 100)## Smoothing the time series...
## Starting wavelet transformation...
## ... and simulations...
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave Phase Ampl Power Power.avg Power.pval Power.avg.pval Ridge Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzwt.image(seg_nee_int.w,
col.contour = "white",
color.palette = mypalette,
lwd = 1,
show.date = TRUE,
main = "US-Seg Flux Tower Net Ecosystem Exchange (gC/m^2)",
plot.ridge = FALSE,
legend.params = list(width = 1.8, label.digits = 4))
Now that we’ve looked at the two variables by themselves, it is time to conduct the cross-wavelet analysis of US-Seg precipitation and NEE:
seg_precip_nee <- seg %>%
select(date, P_int, NEE_int) %>%
as.data.frame()
seg_precip_nee.wc <- analyze.coherency(seg_precip_nee, my.pair = c("P_int", "NEE_int"),
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE,
n.sim = 20)## Smoothing the time series...
## Starting wavelet transformation and coherency computation...
## ... and simulations...## Warning in sqrt(sPower.x * sPower.y): NaNs produced##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave.xy Angle sWave.xy sAngle Power.xy Power.xy.avg Power.xy.pval Power.xy.avg.pval Coherency Coherence Coherence.avg Coherence.pval Coherence.avg.pval Wave.x Wave.y Phase.x Phase.y Ampl.x Ampl.y Power.x Power.y Power.x.avg Power.y.avg Power.x.pval Power.y.pval Power.x.avg.pval Power.y.avg.pval sPower.x sPower.y Ridge.xy Ridge.co Ridge.x Ridge.y Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzwc.image(seg_precip_nee.wc,
main = "Cross-wavelet power spectrum for\nUS-Seg Precipitation and NEE",
legend.params = list(lab = "cross-wavelet power levels"),
periodlab = "period (days)",
col.contour = "white",
color.palette = mypalette,
lwd = .8,
show.date = TRUE,
col.arrow = "red",
which.image = "wp"
)
The cross-wavelet power shows the similarity of the bivariate time series’ wavelet power in the time-frequency domain. Arrows to the right indicate when the two variables are in phase with one another, and arrows to the left indicate when the two variables are in anti-phase.
Final Thoughts
WaveletComp is a useful R package for conducting continuous wavelet-based analyses to identify periodicities in univariate or bivariate time series. Wavelets may be a useful tool for detecting transient high power pulses in a time series. There are many options available in WaveletComp for adjusting a wavelet analysis, and more detailed information can be found in the guided tour and package documentation.
The citation for the WaveletComp R package:
citation("WaveletComp")##
## To cite package 'WaveletComp' in publications use:
##
## Angi Roesch and Harald Schmidbauer (2018). WaveletComp: Computational
## Wavelet Analysis. R package version 1.1.
## https://CRAN.R-project.org/package=WaveletComp
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {WaveletComp: Computational Wavelet Analysis},
## author = {Angi Roesch and Harald Schmidbauer},
## year = {2018},
## note = {R package version 1.1},
## url = {https://CRAN.R-project.org/package=WaveletComp},
## }
##
## ATTENTION: This citation information has been auto-generated from the
## package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.